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Abstract. The Poisson equation occurs in many areas of science and engineering. 
Here we focus on its numerical solution for an equation in d dimensions. In 
particular we present a quantum algorithm and a scalable quantum circuit design 
which approximates the solution of the Poisson equation with error s. The cost is 
almost linear in d and polylog in showing an exponential speedup. The circuit 
uses a number of qubits which is also almost linear in d and polylog in e _1 . We present 
quantum circuit modules together with performance guarantees which can be also used 
for other problems. 



PACS numbers: 03.67.Ac 
1. Introduction 

Quantum computers take advantage of quantum mechanics to solve certain 
computational problems faster than classical computers. Indeed in some cases the 
quantum algorithm is exponentially faster than the best classical algorithm known 
[HElElSlEliniCZlIHlIHl UHl HH Ll2]. It is important to distinguish between two notions 
of "exponentially faster" . We introduce a new dichotomy: 

(i) Exponential quantum speedup: A quantum computer can solve the problem 
exponentially faster than any known classical algorithm. 

(ii) Strong exponential quantum speedup: A quantum computer can solve the problem 
exponentially faster than any classical algorithm. 

The crucial difference between the two concepts is that if a problem satisfies 
criterion 1 someone may invent a faster classical algorithm, while if a problem satisfies 
criterion 2 no faster classical algorithm can exist. 

| Corresponding author. Email: kais@purdue.edu 



Quantum Poisson solver 



2 



To prove strong exponential quantum speedup good bounds on the classical 
computational complexity must be known. Shor's very influential result [8] yields 
exponential speedup but not strong exponential speedup because the classical 
computational complexity of factorization is not known. The same is true of most 
of the quantum speedups for scientific problems. Of course, exponential can be replaced 
by any other function. Thus Grover's algorithm for searching an unstructured database 
[T3] enjoys strong polynomial quantum speedup. 

In this paper we present a quantum algorithm and circuit solving the Poisson 
equation. The Poisson equation plays a fundamental role in numerous areas of science 
and engineering, such as computational fluid dynamics [HI [15], quantum mechanical 
continuum solvation [16], electrostatics [17], the theory of Markov chains [HI EES 120] 
and is important for density functional theory and electronic structure calculations [21 J . 

We will show that for the Poisson equation we achieve strong exponential quantum 
speedup. It is known that the worst case classical complexity of the Poisson equation 
is exponential in the dimension d [22]. The Poisson equation can be solved with 
error e with a number of quantum operations which is almost linear in d and 
polylog in e' 1 . More precisely, the number of quantum operations is proportional 
to max{<i, log 2 e _1 }(log 2 d + log 2 e -1 ) 3 . The quantum circuit uses a number of qubits 
proportional to max{<i, log 2 e _1 }(log 2 d + log 2 e^ 1 ) 2 . 

There are many ways to solve the Poisson equation. We choose to discretize it on 
a regular grid and then solve the resulting system of linear equations. Thus we can use 
the quantum algorithm of Harrow et al. [23J for solving the system without assuming, 
however, that the matrix is given by an oracle. A significant part of our work deals with 
the Hamiltonian simulation of the matrix of the Poisson equation. Moreover, it is an 
open problem to determine whether or not it is possible to simulate a Hamiltonian with 
cost polynomial in the logarithm of the matrix size and the logarithm of £ _1 [21]. Our 
results show that in the case of the Hamiltonian in the Poisson equation the answer is 
positive. 

Our analysis of the implementation includes all the numerical details and will be 
helpful to researchers working on other problems. All calculations are carried out in fixed 
precision arithmetic and we provide accuracy and cost guarantees. We account for the 
qubits, including ancilla qubits, needed for the different operations. We provide quantum 
circuit modules for the approximation of trigonometric functions, which are needed in 
the Hamiltonian simulation of the matrix of the Poisson equation. We show how to 
obtain a quantum circuit computing the reciprocal of the eigenvalues using Newton 
iteration and modular addition and multiplication. We show how to implement quantum 
mechanically the inverse trigonometric function needed for controlled rotations. As we 
indicated, our results are not limited to the solution of the Poisson equation but can 
be used in other quantum algorithms. Our simulation module can be combined with 
splitting methods to simulate the Hamiltonian — A + V, where A is the Laplacian 
and V is a potential function. The trigonometric approximations can be used by 
algorithms dealing with quantum walks. The reciprocal of a real number and a controlled 
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rotation by an angle obtained by an inverse trigonometric approximation are needed for 
implementing the linear systems algorithm [23] regardless of the matrix involved. 

2. Overview 

We consider the (i-dimensional Poisson equation with Dirichlet boundary conditions. 
Definition 1. 

-Au{x) = f{x) x E Id ■= (0, l) d , (1) 
u(x) = x G did, 

where f : Id —> K is a sufficiently smooth function; e.g., see f25\ [2fi{ for details. 

For simplicity we study this equation over the unit cube but a similar analysis 
applies to more general domains in M. d . Often one solves this equation by discretizing it 
and solving the resulting linear system. A finite difference discretization of the Poisson 
equation on a grid with mesh size h, using a 2d + 1 stencil for the Laplacian, yields the 
linear system 

-A h v = f h , (2) 

where fh is the vector obtained by sampling the function / on the interior grid points 
|27[ |26| [28] . The resulting matrix is symmetric positive definite. 

To solve the Poisson equation with error 0(e) both the discretization error and the 
error on the solution of the system should be 0(e). This implies that A/, is a matrix 
of size proportional to e~ ad x e~ ad , where a > is a constant that depends on the 
smoothness of the solution which, in turn, depends on the smoothness of / (291 123 [22] . 
For example, when / has uniformly bounded partial derivatives up to order four then 
a = 1/2. 

There are different ways for solving this system using classical algorithms. Demmel 
[271 Table 6.1] lists a number of possibilities. The conjugate gradient algorithm [30] is 
an example. Its cost for solving this system with error e is proportional to 

e~ ad ^\oge-\ 

where k denotes the condition number of A^. We know k = e~ 2a , independently of d. 
The resulting cost is proportional to e~ ad ~ a loge~ 1 . For details about the solution of 
large linear systems see [31J . Observe that the factor 

£ -ad in 

the cost is the matrix size 

and its contribution cannot be overcome. Any direct or iterative classical algorithm 
solving this system has cost at least e~ ad , since the algorithm must determine all 
unknowns. So any algorithm solving the system has cost exponential in d. In fact 
a much stronger result holds, namely, the cost of any classical algorithm solving the 
Poisson equation in the worst case must be exponential in d [22] . 

We present a scalable quantum circuit for the solution of ([2]) and thereby for the 
solution of Poisson equation with error 0(e) that uses a number of qubits proportional 
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Figure 1: Discretization of the square domain and notation for indexing the nodes. 



to max{d, log 2 e _1 }(log 2 <i + log 2 e -1 ) 2 and a number of quantum operations proportional 
to max{<i, log 2 £ _1 }(log 2 d + log 2 £:~ 1 ) 3 . It can be shown that log 2 d = 0{\og 2 e^ 1 ) 
and the above expressions are simplified to max{<i, log 2 e _1 }(log 2 e' 1 ) 2 qubits and 
max{d, log 2 £~ 1 }(log 2 e^ 1 ) 3 quantum operations. A measurement outcome at the final 
state determines whether the algorithm has succeeded or not. A number of repetitions 
proportional to the square of the condition number yields a success probability 
arbitrarily close to one. 

In Section [3] we deal with the discretization of the Poisson equation showing the 
resulting matrix. We also describe how the matrix in the multidimensional case can 
be expressed in terms of the one dimensional matrix using Kronecker products. This, 
as we'll see, is important in the simulation of the Poisson matrix. In Section H] we 
show the quantum circuit solving the Poisson equation. We perform the error analysis 
and show the quantum circuit modules computing the reciprocal of the eigenvalues and 
from those the controlled rotation needed at the end of the linear systems algorithm 
|23j. In Section [5] we deal with the Hamiltonian simulation of the matrix of the Poisson 
equation. The exponential of the multidimensional Hamiltonian is the d-fold tensor 
product of the exponential of one dimensional Hamiltonian. It is possible to diagonalize 
the one dimensional Hamiltonian using the quantum Fourier transform. Thus it it 
suffices to approximate the eigenvalues in a way leading to the desired accuracy in the 
result. We show the quantum circuit modules performing the eigenvalue approximation 
and derive the overall simulation cost. In Section [6] we derive the total cost for solving 
the the Poisson equation. Section [7] is the conclusion. In Appendix 1 we list a number of 
elementary quantum gates and in Appendix 2 we present a series of results concerning 
the accuracy and the cost of the approximations we use throughout the paper. 

3. Discretization 

3.1. One dimension 

We start with the one- dimensional case to introduce the matrix that we will use 
later in expressing the d- dimensional discretization of the Laplacian, using Kronecker 
products. We have 
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d 2 u(x) 

dx 2 
u(Q) = u{l) 



fix), ie(o,i) 
= o 



(3) 



where / is a given smooth function and u is the solution we want to compute. We 
discretize the problem with mesh size h — 1/M and we compute an approximate solution 
v at M+l grid points Xi — ih, i — 0, . . . , M. Let Ui = u(xi) and fi = f(xi), i = 0, . . . , M. 

Using finite differences at the grid points to approximate the second derivative ([3]) 
becomes 

d 2 u(x) 2ui — Ui^i — u i+ i 



dx 2 



h 2 



6 



where £j is the truncation error and can be shown to be 0(/i 2 ||^f ||oo' 



(4) 

if / has fourth 
derivative uniformly bounded by a constant |27j . 
Ignoring the truncation error, we solve 

/T 2 (-^-i + - v l+1 ) = fi < i < M. (5) 
We have M — 1 equations and M — 1 unknowns f i, t>M-i : 

/ «i \ / 2 -1 \ / v x \ ( h \ 

-1 '•• '•• 



h- 2 -L h 



h 



-2 



\ Vm-1 ) 



\ o 



•. -1 

-1 2 / 



V vm-i ) 



(6) 



V hi-i J 



where Lh is the tridiagonal (M — 1) x (M — 1) matrix above; for the properties of this 
matrix, including its eigenvalues and eigenvectors see (271 Sec. 6.3]. 

3.2. Two dimensions 

In two dimensions the Poisson equation is 



d 2 u(x,y) d 2 u(x,y) 



dx 2 dy 2 
u(x,0) = u(0,y) = 0, x, ye [0,1] 



f(x,y), (x,y)e(0,lf 



(7) 



We discretize this equation using a grid with mesh size h = 1/M; see Figure [TJ 
Each node is indexed Uj t k, j,k G {1,2, ...,M} (Figure [T](a) and (b)). We approximate 
the second derivatives using 

d 2 u u(x — h, y) — 2u(x, y) + u(x + h, y) 

d 2 u , . ^ u(x, y — h) — 2u(x, y) + u(x, y + h) 



dy 2 



h 2 
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Omitting the truncation error, and denoting by — the discretized Laplacian we 
are led to solve 



h 2 + 2v j:k - v j+1)k ) + + 2v j<k - v j)k+1 )) = f j)k 



(8) 



where fj jk = f(jh, kh), j, k — 1, 2, . . . , M — 1 and Vj >k = if j or k e {0, M} i.e., when 
we have a point that belongs to the boundary. 

Using the fact that the solution is zero at the boundary, we reindex ([8]) to obtain 

h~ 2 (Avi - Vi_i - v i+ i - Vi_M+i - Vi+M-i) = fi i — 1,2, ... , (M — l) 2 , (9) 

Equivalently, we denote this system by 

-A h v = f h , 

where is the discretized Laplacian. 

For example, when M = 4, as in Figure dj we have that v = [vi, Vg\ T . 
Furthermore ([9]) becomes 



(10) 





/ v, \ 
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where I is the 3x3 identity matrix, B is 



/ 4 -1 
-1 4 -1 

v - 1 4 



\ 



A is a Hermitian matrix with a particular block structure that is independent of M 
In particular, on a square grid with mesh size h = 1/M we have 

- A h — h~ 2 A 



and A can be expressed in terms of Lh as follows: 

/ L h + 2I -I ••• 
-I L h + 2I -I 

—I 

: '•• '•■ 



A 












V 



(12) 



: : -/ L h + 2I -I 

••• -I L h + 2I J 

and its size is (M - l) 2 x (M - l) 2 [27]. 

Recall that Lh is the (M — 1) x (M — 1) matrix shown in ([6]) and / is the 
(M — 1) x (Af — 1) identity matrix. Moreover, A can be expressed using Kronecker 
products as follows 



A 



I + 1 ®L 



it- 



(13) 
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Figure 2: Overview of the circuit for solving the Poisson equation. Wires with '/' represent 
registers or groups of qubits. W denotes the Walsh- Hadamard transform and 
is applied on every qubit of the register. FT represents the quantum Fourier 
transform. 'HAM-SIM' is the Hamiltonian Simulation subroutine that implements 
the operation e lAt ° . 'INV is the subroutine that computes A -1 . represent 
uncomputation, which is the adjoint of all the operations before the controlled R y 
rotation. 



3.3. d dimensions 



We now consider the problem in d dimensions. Consider the Laplacian 

dxl 



a 

k=i ox k 



We discretize A on a grid with mesh size h—l/M using divided differences. 
As before, this leads to a system of linear equations 

-A h v = f h . (14) 

Note that — = h~ 2 A is symmetric positive definite matrix and A is given by 

d matrices 

and has size (M - l) d x (M - l) d . L h is the (M — 1) x (M - 1) matrix shown in © 
and / is the (M — 1) x (M — 1) identity matrix. See (27] for the details. 
Observe that the matrix exponential has the form 

e M 7 = giL h7 g, e«*7 ( 15 ) 

d matrices 

for all 7 Gffi, where i = \f—l. We will use this fact later in deriving the quantum circuit 
solving the linear system. 



4. Quantum circuit 

We derive a quantum circuit solving the system — A^v = fh, where h = 1/M and 
without loss of generality we assume that M is a power of two. We obtain a solution of 
the system with error 0(e). 
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(i) As in [23] assume the right hand side vector f h has been prepared quantum 
mechanically as a quantum state \fh) and stored in the quantum register B. Nore 
\fh) — Y^j=o ~ l Pj\ u j) where \uj) denote the eigenstates of — A^ and are the 
coefficients. 

(ii) Perform phase estimation using the state \fh) in the bottom register and the unitary 
matrix e - 2m A h /E ^ w h ere j g 2 jt; — [logcf] + log(4M 2 ). The number of qubits in the 
top register of phase estimation is n = 0(}og(E/e)). 

(iii) Compute an approximation of the inverse of the eigenvalues Xj. Store the result 
on a register composed of b = 3|Tog£ _1 ] qubits. The approximation error of the 
reciprocals is at most e. 

(iv) Introduce an ancilla qubit to the system. Apply a controlled rotation on the 
ancilla qubit. The rotation operation is controlled be the register L which stores 
the reciprocals of the eigenvalues of — A/,. The controlled rotation results to 
■sjl — {pd~l Aj) 2 |0) + (Cd/\j)\l), where is a constant. 

(v) Uncompute all other qubits on the system except the qubit introduced on the 
previous item. 

(vi) Measure the ancilla qubit. If the outcome is 1, the bottom register of phase 
estimation collapses to the state Y^jf^o^ _1 PjXj~ 1 \uj) up to a normalization factor, 
where \uj) denote the eigenstates of — A h . This is equal to the normalized solution 
of the system. If the outcome is 0, the algorithm has failed and we have to 
repeat it. An alternative would be to include amplitude amplification to boost the 
success probability. Amplitude amplification has been considered in the literature 
extensively and we do not deal with it here. 

4-1. Error analysis 

We carry out the error analysis to obtain the implementation details. For d = 1 the 
eigenvalues of the second derivative are 

4M 2 sin 2 (j7r/(2M)) j = l,...,M-l. 

For d > 1, the eigenvalues of — Ah are given by sums of the one-dimensional eigenvalues, 
i.e., 

d 

[4M 2 sin 2 (j fc vr/(2M))] 3k = 1, . . . , M - 1, k = 1, . . . , d. 

k=l 

We consider them in non-decreasing order and denote them by Xj, j = 1, . . . , (M — 
l) d . Then Ai = 4dM 2 sin 2 (7r/(2M)) is the minimum eigenvalue and XfM-i) d = 
4dM 2 sin 2 (7r(M - 1)/(2M)) < 4<iM 2 is the maximum eigenvalue. 
Define E by 

\og 2 E= \\og 2 d~\ +log 2 (4M 2 ). (16) 

Then the eigenvalues are bounded from above by E. Recall that we have already 
assumed that M is a power of two. 
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Note that the implementation accuracy of the eigenvalues determines the accuracy 
of the system solution. 

Our algorithm uses approximations Xj, such that |Aj — < < see Theorem 
[2] in Appendix 2. We use n = log 2 E + u bits to represent each eigenvalue, of which 
the log 2 E most significant bits hold each integer part and the remaining bits hold each 
fractional part. Without loss of generality, we can assume that 2 V ^> E. More precisely, 
we consider an approximation A^ of matrix A^ such that the two matrices have the 
same eigenvectors while their eigenvalues differ by at most e. 

We use phase estimation with the unitary matrix e^ lAht °^ E whose eigenvalues are 
e 27TiA,t /(£27r)_ g etting to _ 27r W e obtain the phases (frj = Xj/E G [0, 1). The initial state 
of phase estimation is (Figure 

(M-l) d 

\of n \A) = E &i°n%->> 

where \v,j) is the jth eigenvector of — A h and = (uj\fh), for j = 1, 2, ... , (M — l) d . 
Since we are using finite bit approximations of the eigenvalues, we have 

Xj ^j^ v 
- ~£ - 2 n 

Then <pj 2 n is an integer and phase estimation succeeds with probability 1. 

The state prior to the application of the inverse Fourier transform in phase 
estimation is 

(M-l) d 2 n -l 
j=l k=0 

After the application of the inverse Fourier transform to the first n qubits we obtain 

(M-l) d 

where 

kj = 2 n ( f) j = TXj/E. (18) 

Now we need to compute the reciprocals of the eigenvalues. Observe that 

X t /d = AM 2 sin 2 (vr/(2M)) = 4M 2 (vr/(2M) + 0(M~ 3 )) 2 
= n 2 + 0{M~ 2 ) > 5. 

where the last inequality holds trivially for M sufficiently large. This implies Xj/Cd > 
Xi/C d > 4, where C d = 2^^, for M sufficiently large. We obtain kj = 2 n Xj/E > 
Ai > 4C d . 

Append b qubits initialized to |0) on the left (Reg.L in Figure [2]), to obtain 

(M-l) d 

E /?iion%)h->. 
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Note that from (1181) kj, Xj and Xj/Cd have the same bit representation. The difference 
between the integer kj and the other two numbers is the location of the decimal point; 
it is located after the log 2 E most significant bit in Xj, and after the \og 2 {E/Cd) most 



significant bit in Xj/Cd- Therefore, we can use the labels \kj), 
interchangeably, and write the state above as 

(M-l) d 



Xn ) and 



Xj/C 



Xj/C d )\uj). 



j=i 

Now we need to compute hj := h{Xj/Cd) = Cd/Xj. We do this using Newton iteration. 
We explain the details in Section H~2l We obtain an approximation hj such that 

,2 



< e 



0: 



(19) 



where eo = min{e, -E -1 }. We store this approximation in the register composed of the 
leftmost b = 3 [log 2 eo ] qubits. 
This leads to the state 

(M-l) d 



We append, on the left, a qubit initialized at |0) (Anc. in Figure |2]). We get 

(M-l) d 

/3,|o)|^)|vq)k-). 

0=1 

We need to perform the conditional rotation 

R\0)\u) = (u\l) + Vl -w 2 |0)) \u), 
For this, we will approximate the first qubit by 



< u < 1. 



u'\l) + ^l-(u;') 2 \0), 

with \u — u'\ < el, Ei = min{e, 1/(4M 2 )}. We discuss the cost of implementing this 
approximation in Section 14.31 

The result of approximating the conditional rotation is to obtain 



hj ) , where hj is 



a q = 0(log 2 £: 1 ) bit number less than 1 satisfying \hj — hj\ < e\ and, therefore, 



(20) 



\hj-hj\ <el + el, 

for each j = 1, . . . , (M — l) d . 

Ignoring the ancilla qubits needed for implementing the approximation of the 
conditional rotation, we have the state 

(AI-l) d 



( I — 

3=1 V 



10) 



h. 



Xj/C d )\uj). 
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Uncomputing all the qubits except the leftmost gives the state 

(M-l) d 



E ^ ^ii)+v / i T ^io))io)® 6 iorK.) 



Let Pi = |l)(l| ® I be the projection acting non-trivially on the first qubit. The 

.(j ■ ' 
:'r- 

(M-l) d 



system —AhV = fh has solution Yli=i ^ Pjj'\ u j)- We derive the error as follows 



E ^ii>io> Wn) k->-^i 

j=l 1 



(M-l) d 



a 



(M-l) a 



E ^'pK-) - E 



i=i 

(A/-l) d 



i=i 

(M-l) d 

E E ^ - 1 >. * ■ h j) "j) 



< 




1 1 



1%) 



£ 2 , £ 2 < 1^ + p2 , 2 
fc ' fc l — oi/ ^ fc ' fc l' 



(21) 



where the second from last inequality is obtained using (120]) and the last inequality is 
due to the fact that 
1 1 



A A 



< A-A, A,A>1 



Setting v = [log 2 (17P/£)] gives error £:(l + o(l)) and the number of matrix exponentials 
used by the algorithm is 0(log 2 (E/e)). Therefore, if we measure the first qubit of the 
state and the outcome is 1 the state collapses to a normalized solution of the linear 
system. 



4-2. Computation of X 1 

In this part we deal with the computation of the reciprocals of the eigenvalues, which 
is marked as the TNV module in Figure |2j For this we use Newton iteration to 
approximate v~ l , v > 1. We perform s iterative steps and obtain the approximation x s . 
The input and the output of each iterative step are b bit numbers. All the calculations 
in each step are performed in fixed precision arithmetic. The initial approximation is 
xq = 2~ p , 2 p ~ l < v < 2 P . (We use the notation x% to emphasize that these values have 
been obtained by truncating a quantity Xi to b bits of accuracy, i — 0, . . . , s). 
Theorem [1] of Appendix 2 gives the error of Newton iteration which is 

\x s ~ v' 1 ] <e 2 <e, 

where we have e = min{e , E" 1 } , s = |~log 2 log 2 (2/eQ)] and the number of bits satisfies 
b > 2 \log 2 e^ 1 ] + 0(log 2 log 2 log 2 Eq 1 ). 
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Figure 3: The quantum circuit computing the initial approximation xq 
iteration for approximating v , 2 P ~ 1 < v < 2 P . 



\v) 



2- p of Newton 



Therefore, it suffices that the module of the quantum circuit that computes 1/Xj 
carries each iterative step with 3|"log 2 £o 1 ] qubits of accuracy. 

The quantum circuit computing the initial approximation xq, of the Newton 
iteration is given in Figure [3j The second register holds \v) and is n qubits long, of 
which the first \og 2 {E/Cd) qubits represent the integer part of v and the remaining 
ones its fractional part. The first register is b qubits long. Recall that Xj/Cd > 4. So 
input values below 4 do not correspond to meaningful eigenvalue estimates and we don't 
need to compute their reciprocals altogether; they can be ignored. Hence the circuit 
implements the unitary transformation lO)® 6 ^ ) — > \0)® b \v), if the first \og 2 (E /Cd) — 2 
bits of v are all zero. Otherwise, it implements the initial approximation xq through the 
transformation |0)® 6 |f) — » |xo)|f). 

Each iteration step + 2xi is implemented using a quantum circuit of 

the form shown in Figure H] that computes |£i)|i>) — > \xi+i)\v). This involves quantum 
circuits for addition and multiplication which have been studied in the literature [32] . 

The register holding \v) is n qubits long and the register holding the \xi) and 
\xi + i) is b qubits long. Note that internally the modules performing the iteration 
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— vx, + 2xi 



\v) 



\V) 



Figure 4: Circuit implementing each iterative step of the Newton method. 



steps may use more than b qubits, say, double precision, so that the addition and 
multiplication operations required in the iteration are carried out exactly and then 
return the b most significant qubits of the result. The total number of qubits required 
for the implementation of each of these modules is O(log£o 1 ) an d the total number of 
gates is a low degree polynomial in log^o 1 - 

4-3. Controlled rotation 

We now consider the implementation of the controlled rotation 

R\0)\u) = (u\l) + w 2 |0)) \co), 0<w<1. 

Assume for a moment the we have obtained \9), a q qubit state, corresponding to an 
angle 9 such that sin 9 approximates u. Then we can use controlled rotations R y about 
the y axis to implement R. We consider the binary representation of 9 and have 

q 

d = .e 1 ...e q = Y t e j 2-s, ^e{o,i}. 

Then 

R y {29) = e-»* = ( ^-^ 2d ~ sing 2 ) 
y( \ sin9 x/l-sin 2 ^ J 

= f[e- iYe ^ 3 =f[R e J i (2 1 - J ), 
i=i j=i 

where Y is the Pauli Y operator and 9 e [0, 7r/2]. The detailed circuit is shown in 
Figure 

We now turn to the algorithm that calculates \9) from \ui). Since u corresponds to 
the reciprocal of an approximate eigenvalue of the discretized Laplacian, we know that 
sin -1 (a;) belongs to the first quadrant and sin _1 (u;) = Q(l/M 2 ). Therefore, we can find 
an angle 9 such that | sin(#) — u\ < e\, E\ = min{l/(4M 2 ), e}, using bisection and an 
approximation of the sine function. 

In Appendix 2 we show the error in approximating the sine function using fixed 
precision arithmetic. In Section we show the details of the resulting quantum 
algorithm computing the approximation to the sine function. These results, with a 
minor adjustment in the number of bits needed can be used here. We won't deal with 
the details of the quantum algorithm for the sine function in this section since we present 
them in Section fl5]) that deals with the simulation of Poisson's matrix. We will only 
describe the steps of the algorithm and its cost. 
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|0> 



Ry(l) Ry(l/2) 



Ryil/2"- 1 ) 



\0i) 

|02> 



Figure 5: Circuit for executing the controlled R y rotation. 



Algorithm: 

(i) Take as an initial approximation of 9 the value 7r/4. 

(ii) Approximate the sin(0) with error el/2. Let s e denote this approximation. 

(iii) If s e < to — el/2, set 9 to be the midpoint of the right subinterval. 

(iv) If s e > u + el/2, set 9 to be the midpoint of the left subinterval. 

(v) Repeat the steps 2 to 4 |~log 2 e± 2 ] + 1 times. 

An evaluation at the midpoint of an interval yields a value that satisfies either the 
condition of step 3, or that of step 4, or \sg — u\ < el/2. If at any time both the 
conditions of steps 3 and 4 are false then 9 will not change its value until the end. Then, 
at the end, we have | sin(6>) — u>\ < | sin(6 l ) — sg\ + \sg — uj\ < ef, since the error in 
computing the sine is el/2. On the other hand, if 9 is updated until the very end of the 
algorithm the final value of theta also satisfies | sin(#) — w\ < el, because in the final 
interval we have | sin(0) — u\ < \9 — sin _1 (w)| < e\. 

In a way similar to that of Proposition 1 and Proposition 2 of Appendix 2 we carry 
out the steps of the algorithm in q bit fixed precision arithmetic, q = max{2z/ + 9, 13 + 
v + 2 log 2 M} and sufficiently large v to satisfy the accuracy requirements. (The last 
expression for q is slightly different form that in Proposition 2 because it accounts for 
the fact that in the case we are dealing with here the angle is Vl(l/M 2 )). This gives us 
an approximation to the sine with error 2~ ( - u ~ 1 \ We set 

v = Rog 2 ^ 2 ] + 1. 

Thus v and q are both 0(log 2 £i). 

The algorithm for the sine function is based on an approximation of the exponential 
function using repeated squaring. Each square requires 0(q 2 ) quantum operations and 
0(q) qubits. This is repeated v times before the approximation to the sine is obtained. 
Thus the cost of one bisection step requires 0(z/g 2 ) quantum operations and 0{yq) 
qubits. So, in terms of ei, the total cost of bisection is proportional to (log 2 e 1 -1 ) 4 
quantum operations and (loggej -1 ) 3 qubits. 
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|1> 



b'i) — h 



-ih 2 L h -y 



|1> 



\h) - -A 



-i/i 2 L h 7 



7- 



Figure 6: Quantum circuit for implementing e _lA ' 17 , 76K for the two dimensional discrete 
Poisson equation. The subroutine of e~ lh Lh7 is shown in Figure [71 The registers 
holding \ U2) are m qubits each. 
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?ih~ 2 A-y 



T, 



M 



2M 



T 



M 



-ih 2 L h -y 



Sm-i 

Figure 7: Quantum circuit for implementing e~ th 2Lh "< , 7 e R, where Lh is the matrix in (JS]). 

Sm-i represents sine transform matrix of size (M — 1) x (M — 1), M = 2 m . This 
circuit acts onm + 1 qubits. 



5. Hamiltonian simulation of the Poisson matrix 

In this section we deal with the implementation of the 'HAM-SIM' module (Figure [2]) 
which effectively applies e~ tAht ° onto register B. In our case the eigenvectors of the 
discretized Laplacian are known and we use approximations of the eigenvalues. From 
ffTTj) and (|To| we have 

e -«A h7 = e -ih~ 2 L h y . . . g, e -ih- 2 L h ~f _ , 22 x 
d matrices 

Thus it suffices to implement e~ lh 2L/i7 , for certain 7 e K, 7 = 27r ■ 2 t /E, t = 
0, 1, . . . , log 2 E — 1 that are required in phase estimation. This can be accomplished 
by considering the spectral decomposition SAS of the matrix Lh, where S is the matrix 
of the sine transform [35j [27] . Then S can be implemented using the quantum Fourier 
transform. We will implement an approximation of A. 

We remark that the quantum circuits presented here can be used in the simulation 
of the Hamiltonian — A + V using splitting formulas. For results concerning Hamiltonian 
simulation using splitting formulas see [361 [9]. 

5.1. One- dimensional case 

We start with the implementation of e~ lh 2Lh7 , 7 = 2n2 t /E, E = AM 2 when d = 1 and 
t = 0, 1, . . . , n — 1, where n is the number of qubits in register C; see ffTTj) . The form of 
Lh is shown in ([6]). Lh is a Toeplitz matrix. It is known that this type of matrix can be 
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a) Generic circuit for Tm = Dn, for details refer to [33]. 
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» — x — 



(b) Implementation of B and P m gates in (a). 

Figure 8: Quantum circuit for implementing Tm in Equation [24] and [25] In (b), P m denotes 
the map \x) — > \x-\- \ mod 2 n ) on n qubits. Its implementation is described in 
See Appendix 1 for the definitions of basic gates. 



diagonalized via the sine transform S [37] = SAS', where A is an M — lxM — 1 
diagonal matrix containing the eigenvalues 4 sin 2 (j7r/(2M)), j = 1, . . . M — 1, of Lh and 



S = {Sij}ij=i2,...,M-i is the sine transform where S it j 
Thus 



sin(m*,j=l,---,M-l. 



(23) 



The relationship between the sine and cosine transforms and the Fourier transform can 
be found in [35] Thm. 3.10]. 
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In particular, using the notation in that paper, we have 

Cm+i © (—iSM-i) = 



Cm+i 




(24) 





-iSm-i 

where Cm+i, Sm-\ denote the cosine and sine transforms, and the subscripts M — 1 and 
M + 1 emphasize the size of the respective matrix. F 2 m is the 2M x 2M matrix of the 
Fourier transform. The matrix Tm has size 2M x 2M and is given by (|25p . 



/ 1 



T 



M 



\ 



1 

v/2 



V2 



1 



1 



V2 



V 



1 



(25) 



~V2 



J 



The quantum circuits for implementing the unitary transformation Tm is discussed 
in [33]. The action of T M can be described by [33] 



Tm\0x) = j=\0x) + j^\lx') 
T M \lx) = ^|0x>-^|l^> 



(26) 



where i 2 = — 1, x is an n-bit number ranging 1 < x < 2 n and x' denotes its two's 
complement i.e. x' = 2" — x. The basic idea of implementing Tm is to separate 
its operation into an operator D, which ignores the two's complement in Tm, and a 
controlled permutation tt, which transforms the state \bx) to \bx') only if b is 1. Therefore 
the action of D and 7r can be written as 



(27) 



Clearly, Tm = Dn and the overall circuit for implementing operation Tm is shown in 
Figure 

By ( |24|) the sine transform 5* can be implemented by cascading the quantum circuits 
in Figure El with the circuit for Fourier transform [38 J. An ancilla bit is added to register 
b. It is kept in the state |1) in order to select the block 



D\0x) = 


j= 2 \0x) + 


Til 1 -) 


D\lx) = 


73l 0a; >- 




ti\0x) = 


\0x) 




ir\lx) = 


\lx'} 







-iS M -i 



(28) 
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from the unitary operation T\jF 2 mTm (EI]), a G C. Considering the state \fh), that 
corresponds to the right hand side of (]6]), and for 6, = (i\fh) we have 

(0A,&2,-A/-i) (29) 

values on the 
(M — 1) nodes 

then the element a in Equation [28] has no effect, and the circuit in Figure [6] is equivalent 
to applying {S^^e 2 ™^ /E S M -i) onto the (M - 1) elements of \f h ). This is also 
equivalent to simulating the Hamiltonian e 2mh 2A2t / E with the state 1/^) stored in 
register b. 

We implement e 2mh 2A2t / E where A = {Aj}j=i,...,A/-i is a diagonal matrix 
approximating A = {\j}j=i,... t M-i- 

We obtain each Xj, j = 1, . . . , M — 1 by the following algorithm. 

Eigenvalue Simulation Algorithm (ESA): 

(i) Let r = 2 u+r where v is positive integer which is related to the accuracy of the 
result. The inputs and the outputs of the modules below are s = max{2z/ + 9,11 + 
v + log 2 M} bit numbers. Internally the modules may carry out calculations in 
higher precision 0(s), but the results are returned using s bits. This value of s 
follows from the error estimates in Proposition |2j 

(ii) We perform the transformation 

\j)\0f s ^\j) \y 3 =x J /r) , 

s qubits 

where xj is the s bit truncation of Xj = ji?. Note that yj = Xj/r G (0, 1) and y~j is 
the s bit truncation of yj. Recall that r > 2 and 2M are powers of 2. Calculations 
are to be performed in fixed precision arithmetic, so division does not need to 
be performed actually. All one needs to do is multiply j by n with O(s) bits of 
accuracy, keeping in track the position of the decimal point and then take the s 
most significant bits of the result. 

(iii) We compute the real and imaginary parts of the complex number W\ by truncating, 
if necessary, the respective parts of W = 1 — y 2 + iy to s bits of accuracy; see (H2j) 
in Proposition [TJ This is expressed by the transformation 



Note that since \yj) is s qubits long, Wo can be computed exactly using double 
precision and ancilla qubits and the final result can be returned in s qubits. 
Complex numbers are implemented using two registers, holding the real and 
imaginary parts. Complex arithmetic is performed by computing the real and 
imaginary parts of the result. 



Quantum Poisson solver 

( n - 1 



19 



n — t 
n-t-l 



|0} n -x 
|0)„_ 2 

|0) t _x 
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Figure 9: Quantum circuit for implementing the transformation in Equation 1301 
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(vi) 



We compute W r approximating W{ using repeated squaring. Each step of this 
procedure is accomplished by the transformation 



M(W 2i )) %(W 2J ))\0)® S \0) 



3(W 



2i+ 



which describes the steps in (14211 . The registers holding real and imaginary parts 
of the numbers are s qubits long. 

< ^s(W r ) approximates sin(7rj/(2M)) with error 2~^~ 1 ^ 1 . Hence ^s 2 (W r ) approximates 
the sin 2 (7rj/(2M)). We compute the square of Q{W r ) exactly and multiply it by 
4M 2 (this involves only shifting). We keep the v + log 2 (4M 2 ) most significant bits 
of the result, which we denote by £j. This means that the log 2 (4M 2 ) bits of the 
binary string representing £j compose the integer part and the last v bits compose 
the fractional parts of the approximation to Xj. Then 

For the error estimate details see Proposition [2j When d = 1, n (the number of 
qubits in register C) and v are related by n = v + log 2 (4M 2 ). Moreover, in the one 
dimensional case A, 

For a fixed t, we implement the 



\7 ~~ *-J- 

Let kj be the binary string representing 
transformation 



j ■ 



1%) |o) { 



(30) 



n qubits n qubits 

This is accomplished using CNOTs with the circuit shown in Figure [91 since t <n 
the total number of quantum operations and qubits required to implement the 
circuit for all the values of t is 0(n 2 ). 
(vii) Finally, we use phase kickback (see e.g. [39]) to obtain e 27r *^ 2 ' from the state \kj2 t ) 
where (f)j is the phase corresponding to the eigenvalue £j that approximates Xj] see 

nr 
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5.2. Multidimensional case 

To implement e~ hl , 7 = 2n2 t /E, E denned in ( fT6l) and t = 0, . . . , n — 1 we use 

* . ' V ' 

d matrices 

Therefore the quantum circuit implementing e~ tAhl in d dimensions is obtained by 
the replication and parallel application of the circuit simulating e~ lh 2Lhl . For example, 
when d = 2 we have the circuit in Figure EJ The register B of Figure [2] contains dm 
qubits, m = log 2 M and its initial state is assumed to have the form 

(32) 




values on the nodes of 
(M - l)< x<i ) grid 

where &j = (i\fh)- This way we select SV-i block in tI i F 2 mTm in ( l24"j) in each circuit 
for e~ l?t 2ih7 . Recall that |/^) corresponds to the right hand side of ( |T4"1) . 

The eigenvalues in the d dimensional case are given as sums of the one dimensional 
eigenvalues. We do not need to form the sums explicitly for the simulation of — A^; 
they are computed by the tensor products. The difference between the d dimensional 
and the one dimensional case is that the register C in Figure |2] has |~log 2 d~] additional 
qubits; i.e n = |~log 2 d] + log 2 4M 2 + v. Accordingly, we generate the one dimensional 
approximations to the eigenvalues using the steps 1 — 5 of the eigenvalue estimation 
algorithm of the previous section. Then we append |~log 2 d] qubits initialized to jo)®^ 2 ^ 
to the left of the register holding the and carry out the remaining two steps 6 — 7 
with n = [log 2 d] + log 2 4M 2 + v. The error in the approximate eigenvalues is equal to 
17M 2 d/2 u ; see Theorem 



5.3. Simulation cost 

Simulating the sine and cosine transforms ff24|) requires 0(m 2 ), m = log 2 M quantum 
operations and 0(m) qubits [33]. The diagonal eigenvalue matrix of the one dimensional 
case ( J23|) is simulated by ESA. Its steps 1—3 and step 5 require 0(s 2 ) quantum operations 
and O(s) qubits. In step 4 repeated squaring is performed v + 7 times. Each repetition 
or step of the procedure requires 0(s 2 ) quantum operations and O(s) qubits. The total 
cost of step 4 is proportional to v ■ 0(s 2 ) quantum operations and v ■ O(s) qubits, 
accounting for any ancilla qubits used in repeated squaring. Step 6 requires 0{n + 1) 
quantum operations and qubits for fixed t. Step 7 requires 0(n 2 ) quantum operations, 
due to the Fourier transform, and 0{n) qubits. 

Using Theorem [2} and requiring error e in the approximation of the eigenvalues, we 
have So we set 

17 E 



v = Rog 2 — 1 , 
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i.e. v = B(log 2 d + m + log 2 e 1 ). We also have n = Q{v) and s = G(n). 
We derive the simulation cost taking the following facts into account: 

• Steps 1 — 5 deal with the approximation of the eigenvalues. These computations 
are not repeated for every t — 0, . . . , n — 1. The total cost of these steps is 0(n 3 ) 
quantum operations and 0(n 2 ) qubits. 

• The total cost of step 6, resulting from all the values of t, is 0(n 2 ) quantum 
operations and qubits. 

• The total cost of step 7, that applies phase kickback for all values of t, does not 
exceed 0{n 3 ) quantum operations and 0{n 2 ) qubits. 

Therefore the total cost to simulate e~ lh Lhl , 7 = 2n2 t /E, for all t — 0, . . . , n — 1, is 
0{n 3 ) quantum operations and 0{n 2 ) qubits. From ( 1221) we conclude that the cost to 
simulate Poisson's matrix for the d dimensional problem is d-0(n 3 ) quantum operations 
and d ■ 0{n 2 ) qubits. 

Finally, we remark that the dominant component of the cost is the one resulting 
from the approximation of the eigenvalues (i.e., the cost of steps 1 — 5). 

6. Total cost 

We now consider the total cost for solving the Poisson equation (CQ). Discretizing the 
second derivative operator on a grid with mesh size h = 1/M results to a system of 
linear equations, where the coefficient matrix is (M — l) d x (M — l) d , i.e. exponential in 
the dimension d > 1. Solving this system using classical algorithms has cost that grows 
at least as fast as the number of unknowns (M — l) d . For the case d = 2, [27] Table 6.1] 
summarizes the cost of direct and iterative classical algorithms solving this system. 

For simulating Poisson's matrix we need d-0(n 3 ) quantum operations and d-0(n 2 ) 
qubits, where n = 0(log 2 d + m + log 2 e -1 ) and m = log 2 M. To this we add the cost for 
computing the reciprocal of the eigenvalues which is 0((log 2 Bq 1 ) 2 log 2 log 2 Eq 1 ) quantum 
operations and 0((log 2 e$ x ) log 2 log 2 Sq x ) qubits, accounting for the 0(log 2 log 2 Sq x ) 
Newton steps, eo = mm{e, E' 1 }. Finally, we add the cost of the conditional rotation 
which is proportional to (log 2 e^ 1 ) 4 quantum operations and (log 2 e^ 1 ) 3 qubits, E\ = 
min{l/(4M)V}. 

From the above we conclude that the quantum circuit implementing the algorithm 
requires of order d ■ 0(n 3 ) + (log 2 £:^ 1 ) 4 quantum operations and d ■ 0(n 2 ) + (log 2 e|f 1 ) 3 
qubits. 

The relation between the matrix size and the accuracy is very important in assessing 
the performance of the quantum algorithm solving a linear system, since its cost depends 
in both of these quantities [23J . In particular, for the Poisson equation we have ignored, 
so far, the effect of the discretization error of the Laplacian A. If the grid is too coarse 
the discretization error will exceed the desired accuracy. If the grid is too fine, the 
matrix will be unnecessarily large. Thus the mesh size and, therefore, the matrix size 
should depend on e, i.e. M = M(e). This dependence is determined by the smoothness 
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of the solution u, which, in turn, depends on the smoothness of the right hand side 
function /. For example, if / has uniformly bounded partial derivatives up to order 
four, then the discretization error is 0(h 2 ) and we set M = e^ 1 ^ 2 ; see [271 [26] f° r details. 
In general, we have M = e~ a , where a > is a parameter depending on the smoothness 
of the solution. This yields n = 0(log 2 <i + log 2 e _1 ), since m = log 2 M = alog 2 £ _1 . 
The resulting number of the quantum operations for the circuit is proportional to 

max{d, log 2 e _1 }(log 2 d + log 2 e -1 ) 3 , 

and the number of qubits is proportional to 

max{d, log 2 e _1 }(log 2 rf + log^ 1 ) 2 , 

It can be shown that log 2 G? = 0(log 2 £ _1 ) and the number of quantum operations and 
qubits become proportional to 

max{d, log 2 e~ 1 } (log 2 e~ 1 ) 3 , 

and 

max{d,log 2 e _1 }(log 2 £: _1 ) 2 , 

respectively. 

Observe that the condition number of the matrix is proportional to e~ 2a and is 
independent of d. Therefore a number of repetitions proportional to e~ Aa leads to a 
success probability arbitrarily close to one, regardless of the value of d. In contrast 
to this, the cost of any deterministic classical algorithm solving the Poisson equation 
is exponential in d. Indeed, for error e the cost is bounded from below by a quantity 
proportional to e~ d l r where r is a smoothness parameter [22J. 

7. Conclusion 

We present a quantum algorithm and a circuit for approximating the solution of the 
Poisson equation in d dimensions. The algorithm enjoys strong exponential quantum 
speedup and the circuit is scalable. We also provide quantum circuit modules and 
performance guarantees which can be used for other problems. The modules include 
reciprocal of eigenvalues and trigonometric approximations. 

Acknowledgements 

We would like to thank NSF CCI center, " Quantum Information for Quantum Chemistry 
(QIQC)", Award number CHE-1037992 and Army Research Office (ARO) for financial 
support. 



Quantum Poisson solver 



23 



Appendix 1 

In this paper, X, Y and Z are Pauli matrices o~ x , a y and a z . I represents identity matrix. 
H is the Hadamard gate and W in Figure [2] represents H® n where n is the number of 
qubits in the register. The matrix representations of other quantum gates used are the 
following: 





RzM = e i0 [ ; " , R z {6) 





„i6 



RM ~ 1 zsin(f) cos(f) ) ' ^ " [ sin(f) cos(|) ' (34) 




Appendix 2 

Theorem 1. Consider the approximation x s to v~ l , v > 1, using s steps of Newton 
iteration, with initial approximation x = 2~ p , 2 P_1 < v < 2 P . Assume that each step 
takes as inputs b bit numbers and produces b bit outputs and that all internal calculations 
are carried out in fixed precision arithmetic. Then the error is 

\x s — v _1 | < En + s2~ b , 

where En denotes the desired error of Newton iteration without considering the 
truncation error, En > 2~ 2 ° . The truncation error is given by the second term and 
s > [log 2 log 2 e^ 1 ] , b> p. 

Proof. Consider the function g(x) — 1/x — v, x > 0, where g(l/v) = 0. The Newton 
iteration for approximating the zero of g is given by 

x s+1 = Lp(x s ) = 2x s -vx] s = 0, 1, . . . . 

The error e s = \x s — l/v\ satisfies e s+ \ = ve 2 s . Unfolding the recurrence we get 

e s < (ve ) 23 . 

Now consider the least power of two that is greater than or equal to v, i.e., 2 P ~ 1 < v < 2 P . 
Clearly p > 1 since v > 1 and veo < 1/2. For error en we have 2~ 2S < en, which implies 

S > [log 2 log 2 Ejf 1 ] . 

The derivative of the iteration function is decreasing and we have |<//| < 2(1— axo) < 
1. We will implement the iteration using fixed precision arithmetic. We first calculate 
the round off error. We have 



x = x 

&i = (p(x ) + ei 
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x s = </?(x s _i) + e s , 
where the ej denotes truncation error at the respective steps. Thus 

x s - x s = <^(x s _i) + e s - (p(x s -i), 
and using the fact \ip'\ < 1 we obtain 

s 

\x s X s | ^ |^s— 1 %s~ l| |Cs| ^ ^ ^ ^ s2 , 

i=l 

assuming that we truncate the intermediate results to b bits of accuracy. 



□ 



Lemma 1. Let x G [tt/(2M), tt/2) and W = 1 + if - 



%. T/ien 



/or r > 2. 
Proof. e lx 



W r \ < r~\ 



(e lx l r ) r = (W + E(x/r)) r , where for y = x/r, E(y) = 



£ 

fc>3 



fc! 



< 



fc-3 



fc>3 



fc>3 



fc! 



x 



£ 

fc>0 



Ml 

fc>3 fc! 

fc! 



and 



1 2/1 



(A; + 3)! fc! 
^^<x (36) 

(fc+l)(fc+2)(fc+3) - 6 



< M! e ly| < 13,13 
— 6 11 

where the last inequality holds for \y\ = |-| < 1, which is true due to our assumptions. 
Hence < |f | 3 for \x\ < r. 

We then turn our attention to the powers of W. 



2 2 

i _ _ _ i . JU JU . tAj tXj 

W = 1 + i ? < 1 + - + — 



(37) 



For all fc G {l,2,...,r} we have, 
|^| fc <(! + - + 



1 X !V< e (^) fe 
r r 2 / 



e^e^r* < e N e ^- < e 2. < e * (3 8 ) 



where we have used the fact that * < 1. The second inequality is due to (1 + a) fc < e ka , 
a G R, fc G Z+ Indeed 



(! + «)* = £ 

fc 

£ 



2=0 



.fc-Z 



£ 

2=0 



fc! 



' = £ 



fc! (fc a y 



Z!(fc-Z)! ^ «(*-*)! 



k k ~ l 



k fc(fc - 1) • • • (/ + 1) / • • • 1 (fca)*"' < ^ (fca 



fc-/ 



«=o v - 



fc fc_z 



(39) 



< i 



< i 



Quantum Poisson solver 

Finally we look at the approximation error. Note that 



25 



r 
k 



W r + 



W r ~"E [-) + 



W k [E (f )] 
W° 



i — k 



E\ X - 

r 



(40) 



error in r-th power 

Consider the fc-th term in the error series. According to (1561) we have 



r 
k 



\W\ r 



r—k 



E[- 

r 



< c 

= C 

< C 



r \ .x 
k } r 



3k 



C 



v ' 



\x 



3k 



k\(r — k)\ r 
r(r — 1) • • • (r — k + 1) 1 \x 



?,k 



3k 



k\ r 2k 



7T 



2fc . , 

7T „ f W 

< -eM — 



2k 



where C = e 71 and we use Stirling's formula k\ = \/2~rk k+1 ^ 2 exp f— A; + jf-), # £ (0, 1), 
[301 p. 257] to obtain |a;| fc /A;! < h~ k x k e k < 1 for A; > 5, since |x| < |. So the total 
approximation error is bounded by 



\ e ix -W r \ <^{r)k\W 



r—k 



k=l 



3k 77- 

~ 2 v 



\x\ 
r 



-<2 7 --(41) 
□ 



Lemma 2. Under the assumptions of LemmaU\ 
\smx-%(W r )\ < 2 7 /r 

and 

\cosx-^(W r )\ < 2 7 /r. 

The proof is trivial and we omit it. 

Proposition 1. Let r = 2 V+1 for v > 1 and consider the procedure computing W r , 
as defined in Lemma [I] using repeated squaring. Assume each step computing a square 
carries out the calculation using fixed precision arithmetic and that its inputs and outputs 
are s bit numbers. Let W r be the final result. Then the error is 

2^+9 



W r -W r 



< 



for s > ll + z/ + log 2 M, where 1/M is the mesh size in the discretization of the Poisson 
equation. 
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Proof. We are interested in estimating sin(j7r/(2M)), for j = 1,2, . . . , M — 1. We 
consider x G [it/(2M),tt/2). We approximate e lx and from this sin a;; the imaginary 
part of e lx . Let y = f- < 2~ 7 . We truncate it to s bits of accuracy to obtain y. Note that 
W = 1 - y 2 + iy satisfies \W\ 2 = 1 - y 2 + y 4 < 1. Let W = 1 - y 2 + iy, \y - y\ < 2~ s . 
Then (Wol 2 < \W\ 2 + Ay2~ s < 1, for s > 11 + v + \og 2 M. This value of s follows by 
solving 

Ay2- S < y 2 /2, 
which ensures that Wq < 1. In addition 



Sft Wn - W 



and 



<2y2~ 



< 2~ 



-2* 



Define the sequence of approximations 
Wi = W + ei 



Wo 



W 2 + e 2 



W r 



r/2 



+ e r 



(42) 



where r = 2 U+7 and the error terms e±, e%, . . . , e r are complex numbers denoting that 
the real and imaginary parts of the results are truncated to s bits of accuracy. 

Observe that if IW 7 ^- 1 ! < 1 then \W 2 j\ < 1, since IW^j- 1 ! 2 < 1 and truncation of 
real and imaginary parts does not increase the magnitude of a complex number. Since 
| Wo I < 1, all the numbers in the sequence (j4"2j) belong to the unit disk S in the complex 
plane. 

Let z = a + hi. Then the function that computes z 2 can be understood as a vector 
valued function of 2 variables, h : 5* — )■ S, such that h(a,b) = (a 2 — b 2 ,2ab). The 
Jacobian of h is 

a —b 
b a 



J 



(a, b) e S 



and its Euclidean norm satisfies || J|| < 2, since a 2 + b 2 < 1. Using this bound we obtain 

\W r - W r \ < \W r - (W r/2 ) 2 \ + |e r | 

< 2{2\W r/4 -W r/4 \ + |e r/4 |} + \e r \ 

< 2 V+7 \W - Wi| + 2 u+7 - 1 \e 2 \ + ... + 2°\e 2 ,+7\ 



■tv+7 



< 2 



u+7 



W -Wo 
W -W 



+ 2' /+7 \e 1 \ + . . . + \e 2 »+7\ 



+ 



V2 



v+7 
j=0 



+7-3 
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< 2 



< 4- 



H-7, 



1 y 1 

22*/ + 2^ 



V2 



(2^+ 8 - 1) 



where the last inequality follows since 2y + 2 s < 2 6 + 2 



-ii 



(43) 
□ 



Proposition 2. Under the assumptions of PropositionUlwe approximate smx by ^s(W r ), 
x E [n/(2M),n/2), with s = max{2z/ + 9, 11 + u + log 2 M} bits and r = 2 V+1 . Then the 
error is 



sin 



x-S(Wr)\ < 2- { - v - l) . 



Moreover: 



Denoting by W r j the approximations to sin (717/ (2M)) ; j = 1, 2, . . . , M — 1, we have 
the following error bound 

AM 2 sin 2 (j7r/(2M)) - AM 2 (^{W r ,j 



< 2-^M 2 , 



j = 1, 2, . . . , M — 1, for the eigenvalues of the matrix h 2 Lh that approximates the 
second derivative operator, using mesh size h = 1/M. 

Letting £j be the truncation of AM 2 ( Q(W r j) ) to v bits after the decimal point (the 



length of lj is v + log 2 (4M ) bits, and v is sufficiently large to satisfy the accuracy 
requirements) we have 



\AM 2 sin 2 (j7r/(2M)) - tA < 17 • 2""M 2 



forj = l,2,...,M-l. 
Proof. We have 

e ix -W r < \e ix -W r \+ W r -W r 

2 7 2 U+9 

< 1 

- 2^+7 2 s 

= 2~" + = 7, (44) 

2 s 2 V 

for s = max{2z/ + 9, 11 + v + log 2 M}, which completes the proof of the first part. The 
proof of the second and third part follows immediately. □ 



Theorem 2. Consider the eigenvalues 

d 



A 



31,— ,3d 



AM 2 Y[ si 



sin 



fc=i 



JfcTT 

2M 



j k = 1, 2, . . . , M - 1, k = 1, 2, . . . , d of - A h , h = 1/M. Let 



A 



3l,— ,3d 



k=i 
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where £j k are defined in Proposition^ ji~ = 1, 2, . . . , M — 1, k = 1, 2, . . . , d. Then 

The proof follows from Proposition [2] and the fact that the d dimensional eigenval- 
ues are sums of the one dimensional eigenvalues. 
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